library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3
library(neonAOP)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(rhdf5)
library(rgeos)
## rgeos version: 0.3-19, (SVN revision 524)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-3
## Polygon checking: TRUE
library(ggplot2)
chm <- raster("../NEONdata/D17-California/SOAP/2013/lidar/SOAP_lidarCHM.tif")
We wanted to match up our in situ data with the hyperspectral data so that we could subset a hyperspectral flight line. As a first step, we identified plot boundaries based on the stem locations in the in situ data.
stem.map <- readOGR("../NEONdata/D17-California/SOAP/2013/insitu/veg-structure",
"soap_stems")
## OGR data source with driver: ESRI Shapefile
## Source: "../NEONdata/D17-California/SOAP/2013/insitu/veg-structure", layer: "soap_stems"
## with 1231 features
## It has 26 fields
# look at a plot
plot(chm)
plot(stem.map, add=TRUE)
# group stems by plotid, record the max and min northing and easting values
# this will be used later to create a shapefile for plot boundaries
stem.map.extent <- stem.map@data %>%
group_by(plotid) %>%
summarise(northing.max = max(northing) + 5,
northing.min = min(northing) - 5,
easting.max = max(easting) + 5,
easting.min = min(easting) - 5)
# assign new variables for use with previously created code
yPlus <- stem.map.extent$northing.max
yMinus <- stem.map.extent$northing.min
xPlus <- stem.map.extent$easting.max
xMinus <- stem.map.extent$easting.min
# code from NEON tutorial on creating square plot extents
square <- cbind(xMinus, yPlus,
xPlus, yPlus,
xPlus, yMinus,
xMinus, yMinus,
xMinus, yPlus)
ID <- stem.map.extent$plotid
# Create a function to do this
polys <- SpatialPolygons(mapply(function(poly, id) {
xy <- matrix(poly, ncol=2, byrow=TRUE) # take a list and create a matrix
Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID),proj4string=CRS(as.character("+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")))
polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))
plot(chm)
plot(polys.df, add=TRUE)
We used Leah’s code to look through all the hyperspectral flightlines. I’m not going to rehash it all here, but we decided to narrow down to one flightline that covered 4 plots.
We also used Leah’s code to get extents for all flightlines and saved this on our local computer.
## SOAP Clip
# the name of the site
site <- "SOAP"
domain <- "D17"
fullDomain <- "D17-California"
level <- "L1"
dataType <- "Spectrometer"
level <- paste0(site,"_L1")
year <- "2013"
productType <- paste0(site,"_", dataType)
dataProduct <- "Reflectance"
drivePath <- "Volumes"
driveName <- "AOP-NEON1-4"
dataDir <- file.path(drivePath, driveName,
domain,
site, year, level, productType, dataProduct)
dataDir <- paste0("/", dataDir)
The right boundary of the flightline appears on the plot below
flight1 <- readOGR("exports/SOAP_flightLines","NIS1_20130612_104651_atmcor")
## OGR data source with driver: ESRI Shapefile
## Source: "exports/SOAP_flightLines", layer: "NIS1_20130612_104651_atmcor"
## with 1 features
## It has 1 fields
# look at this with our plots
plot(chm)
plot(polys.df, add=TRUE)
plot(flight1, add=TRUE)
Note that 4 plots are close to the center of the flightline. We will subset for those 4 plots.
# choose the plots that intersect with flight 1 for extracting HSI
flight1.plots <- intersect(polys.df, flight1)
flight1.plots
## class : SpatialPolygonsDataFrame
## features : 4
## extent : 297029, 297132.1, 4100078, 4100919 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 2
## names : id.1, id.2
## min values : SOAP139, NIS1_20130612_104651_atmcor.h5
## max values : SOAP555, NIS1_20130612_104651_atmcor.h5
# check this subset
plot(chm)
plot(flight1.plots, add=TRUE)
Now we create a boundary that includes all 4 plots to subset the hyperspectral data.
# thanks for the code, leah!
# define the CRS definition by EPSG code
epsg <- 32611
# define the file you want to work with
# this is the hyperspectral flightline from the hard drive
f <- paste0(dataDir, "/NIS1_20130612_104651_atmcor.h5")
# define clip.extents
clip.extent <- flight1.plots
# calculate extent of H5 file
h5.ext <- create_extent(f)
h5.ext
## class : Extent
## xmin : 296670
## xmax : 297409
## ymin : 4098356
## ymax : 4103734
# turn the H5 extent into a polygon to check overlap
h5.ext.poly <- as(extent(h5.ext),
"SpatialPolygons")
crs(h5.ext.poly) <- crs(clip.extent)
# test to see that the extents overlap
gIntersects(h5.ext.poly,
clip.extent)
## [1] TRUE
# Use the clip extent to create the index extent that can be used to slice out data from the
# H5 file
# xmin.index, xmax.index, ymin.index, ymax.index
# all units will be rounded which means the pixel must occupy a majority (.5 or greater)
# within the clipping extent
index.bounds <- vector("list", length(clip.extent))
index.bounds <- calculate_index_extent(extent(clip.extent),
h5.ext)
index.bounds
## [1] 359 462 2815 3656
# this is what i wrote to a csv!